Sexual dimorphism in the tardigrade Paramacrobiotus metropolitanus transcriptome

Background In gonochoristic animals, the sex determination pathway induces different morphological and behavioral features that can be observed between sexes, a condition known as sexual dimorphism. While many components of this sex differentiation cascade show high levels of diversity, factors such as the Doublesex-Mab-3-Related Transcription factor (DMRT) are widely conserved across animal taxa. Species of the phylum Tardigrada exhibit remarkable diversity in morphology and behavior between sexes, suggesting a pathway regulating this dimorphism. Despite the wealth of genomic and zoological knowledge accumulated in recent studies, the sexual differences in tardigrades genomes have not been identified. In the present study, we focused on the gonochoristic species Paramacrobiotus metropolitanus and employed omics analyses to unravel the molecular basis of sexual dimorphism. Results Transcriptome analysis between sex-identified specimens revealed numerous differentially expressed genes, of which approximately 2,000 male-biased genes were focused on 29 non-male-specific genomic loci. From these regions, we identified two Macrobiotidae family specific DMRT paralogs, which were significantly upregulated in males and lacked sex specific splicing variants. Furthermore, phylogenetic analysis indicated all tardigrade genomes lack the doublesex ortholog, suggesting doublesex emerged after the divergence of Tardigrada. In contrast to sex-specific expression, no evidence of genomic differences between the sexes was found. We also identified several anhydrobiosis genes that exhibit sex-biased expression, suggesting a possible mechanism for protection of sex-specific tissues against extreme stress. Conclusions This study provides a comprehensive analysis for analyzing the genetic differences between sexes in tardigrades. The existence of male-biased, but not male-specific, genomic loci and identification of the family specific male-biased DMRT subfamily provides the foundation for understanding the sex determination cascade. In addition, sex-biased expression of several tardigrade-specific genes which are involved their stress tolerance suggests a potential role in protecting sex-specific tissue and gametes. Supplementary Information The online version contains supplementary material available at 10.1186/s40851-024-00233-0.


Introduction
Reproductive modes in animals are typically categorized into two major categories: asexual and sexual.Sexually reproducing animals produce sex-specific gametes, and genetic exchange between sexes leads to higher genetic diversity [1,2].Gonochoristic animals usually show sexual dimorphism, not only in gametes, but also in somatic tissues, physiology, and behavior within a single species, demonstrating dynamic intraspecies differentiation.
Although many aspects of the mechanisms for inducing sex differences remain to be elucidated, they are usually regulated by sex-specific differences in the genome and gene expression [3].Gonochoristic animals must undergo sex determination through common wellstudied mechanisms that underlie the development of sex-specific organs.The physiological systems of sex determination vary among species, but are generally categorized into two types: determination by sex-linked chromosomes and environmental cues [4,5].Species possess sex chromosomes that show different karyotypes depending on their sex.In contrast, environmental cues, including temperature, nutritional status, and population density, act as initial cues for sex determination [6].Regardless of the mode of sex determination, several widely conserved genes play crucial roles in sex-specific organ development.
The transcription factor family Doublesex and Mab-3 Related Transcription factor (DMRT) comprises key regulators of somatic tissue development in various animals [7].In animals utilizing sex chromosomal sex determination systems, not only DMRT orthologs on the sex chromosome, but also on the autosomes, are involved in sex determination cascades to regulate the growth of sex-specific tissues [8,9].In contrast to the chromosomal sex determination system, environmental cues induce the development of sex-specific tissues in normally parthenogenetic individuals through the expression of DMRT orthologs (e.g., Dsx1 in Daphnia magna), leading to genetic exchange through mating [10][11][12].The interplay between highly diverse and conserved components in generating different sexes to overcome environmental and genetic challenges presents a significant challenge to the understanding of the sex determination cascade.
The phylum Tardigrada, a member of the Ecdysozoa with over 1,400 species [13], is divided into three classes: Heterotardigrada, Eutardigrada, and nomen dubium Mesotardigrada [14,15].Tardigrades are renowned for their ability to tolerate extreme environments, and studies have identified tardigrade-specific proteins that mediate tolerance against nearly complete desiccation and anhydrobiosis using species from the family Echiniscidae (Heterotardigrada) [16], Hypsibiidae (Eutardigrada) [17], and Macrobiotidae (Eutardigrada) [18].Asexual (parthenogenesis) and sexual reproduction have been observed within this phylum, with reported instances of both gonochorism and hermaphroditism in sexually reproducing species [19].Sexual dimorphism in morphology and behavior during mating have been observed [20,21].In contrast, we lack knowledge on the molecular mechanisms that induce sexual dimorphism, as most molecular and genomic studies have focused on parthenogenetic species [22].
To this end, we conducted genomic and transcriptomic comparisons between males and females of the model gonochoristic tardigrade Paramacrobiotus metropolitanus to identify the molecular factors related to sexual dimorphism.This species, which is rich in ecological information, has a reported 170 Mbp genome, is relatively easy to culture, and shows a male-biased sex ratio (Male:Female = 7:3), but morphological sexual dimorphism, excluding testis/ovary, has not been described [18,21,[23][24][25].The results in this study lay the foundation for subsequent studies aimed at identifying a master regulator of the sex determination cascade and sex-dependent genetic differences in tardigrades.

Tardigrade culture condition and specimen preparation
The tardigrade P. metropolitanus TYO strain was cultured following methods described in a previous report [24].The specimens were sexed using the method described by Sugiura et al. [23].The eggs of P. metropolitanus were individually placed in an agar-coated dish, and hatched individuals were separated and reared separately to avoid sex contamination.These specimens were then grown until the development of sexual organs that were used for sexing.

RNA extraction and sequencing
Total RNA was extracted as described by Arakawa et al. [26].Two hundred and fifty specimens of each sex were placed in a 1.5 mL tube with minimal water, and 100 µL of TRIzol reagent was added (Thermo Fisher Scientific).Total mRNA was extracted using the Direct-zol RNA kit (Zymo Research), and the samples were transported to Chemical Dojin for sequencing.The transcriptome sequencing libraries were prepared with poly A selection using the NEBNext Ultra II RNA Library Prep Kit for Illumina (New England Biolabs) and were sequenced using the NovaSeq 6000 instrument (Illumina, 150 bp PE).Four and three replicates were prepared for males and females, respectively.
To extract genomic regions enriched in transcripts biased to either sex, we performed enrichment analysis based on the number of DEGs with more than 10 × fold change within 200 kbp windows (100 kbp steps).Genomic bins were created with BEDtools v2.31.1 [43] and the number of genes fitting the criteria was calculated using BEDtools intersect.An in-house Rscript was used to perform Fisher's exact test for each bin, and p-values were corrected by BH method.Regions with Q-value < 0.01 were considered enriched.

Phylogenetic analysis
To identify and analyze the expression patterns of DMRT genes, we first performed an exhaustive search for genes harboring Doublesex-Mab-3 related domains (DM domains).Initial candidates were extracted based on the InterProScan searches performed above, and the corresponding amino acid sequences were submitted to a BLASTP v2.2.22 [49] search against P. metropolitanus.In addition, the amino acid sequences of P. metropolitanus DMRT orthologs were subjected to a BLASTP search against amino acid sequences predicted from various tardigrade genomes [33].The amino acid sequences for the tardigrade DMRT orthologs, metazoan orthologs provided in a previous study [50], and velvet worm DMRT ortholog were pooled (Table S1) and then aligned with MAFFT v7.450 [51] and subjected to phylogenetic tree construction using IQTREE2 v2.2.2.6 [52].The phylogenetic tree was visualized in FigTree v.1.4.3 (http:// tree.bio.ed.ac.uk/ softw are/ figtr ee).The expression patterns of H. exemplaris and R. varieornatus DMRT orthologs during developmental stages were obtained from a previous report by our group [53].Additional alignments for the DMRT proteins were performed by MAFFT v7.450 and visualized using MView (https:// www.ebi.ac.uk/ Tools/ msa/ mview/).
Similarly, we conducted a phylogenetic analysis of CAHS genes.We obtained annotated CAHS sequences from our previous report [54] and pooled the amino acid sequences of P. metropolitanus CAHS candidate orthologs [18].A phylogenetic tree was constructed using the same procedure.

Genome extraction
Virgin P. metropolitanus was prepared by the method described above, and a single tardigrade was placed in a 0.2 ml tube after 1% penicillin/streptomycin treatment for 2 h to remove contamination.Genomic DNA was extracted and prepared using the method described by Arakawa et al. [26].An individual was crushed by pressing it against a tube wall using a pipette tip.Genomic DNA was extracted with Quick-gDNA MicroPrep kit (Zymo Research) with three freeze-thaw cycles and then following the manufacturer's protocol.The extracted DNA was sheared to 550 bp target fragments with Covaris M220 and an Illumina library was constructed with a Thruplex DNA-Seq kit (Takara BioRubicon Genomics).Quantification, quality, and library size selection were performed with Qubit Fluorometer (Life Technologies) and TapeStation D1000 ScreenTape (Agilent Technologies), respectively.Sequencing library fragments in the range of 400-1,000 bp were cut and purified with a NucleoSpin Gel and PCR Clean-up kit (Clontech) and sequenced using a NextSeq500 sequencer with HighOut-putMode 75 cycles kit (Illumina).The reads were de-multiplexed, and adaptor sequences were removed using the bcl2fastq v2 software (Illumina).

Sex specific region analysis
To identify candidate sex chromosome regions, we employed the Y chromosome genome scan (YGS) method [68], which was previously used to identify Drosophila melanogaster sex chromosome contigs.Briefly, reads from the same sex were pooled, and 15-mers were extracted with jellyfish count v2.2.4 or v2.2.10 [69].Scripts from the YGS method v.11b (8 Oct 2012 10AM) were then used to calculate the percentage of validated single-copy unique k-mers (P_VSC_UK) for each contig.This was performed for the previously published genome, as well as for the SPADES assembly performed above.We also tested the coverage for both sexes calculated from the gDNA-Seq data.The raw gDNA-Seq reads were mapped to the genome using BWA-mem2 v2.2.1 and converted into BAM files using SAMtools v1.17.The genome was split into 10 kbp bins and the average coverage for each bin was calculated using BEDtools v2.31.0.The values were then normalized by the median of all bins for that sample, and the average for males and females was computed.

Genotyping for male specific regions
Virgin specimens were replaced with single-worm lysis buffer (50 mM KCl, 10 mM Tris pH 8.2, 2.5 mM MgCl2, 0.45% NP-40, 0.45% Tween20, 0.01% gelatin, 2 μg of Proteinase K) [73].The specimen was then dissolved by freeze-thaw cycles (three times for liquid N 2 and RT) and incubated at 60 °C for 1.5 h and 95 °C for 25 min.Genotyping PCR was performed using the following conditions: 94 °C for 3 min; 40 cycles of 94 °C for 30 s, 50 °C for 30 s, and 68 °C for 1 min; and a final extension at 68 °C for 5 min.Primer sequences were designed using Primer3 [74] from the nucleotide sequences of scaffold Parri_scaffold0000295 (Table S2).Quick-Taq (TOYOBO) was used for the polymerase with the concentrations of each reagent, following the manufacturer's instructions.Electrophoresis was performed at 100 V for 20 min with 1.2% agarose gel/TAE (NacalaiTesque), and then the gel was strained with ethidium bromide for 20 min.The DNA bands were visualized using ChemiDoc (BioRad).

Transcriptomics of P. metropolitanus sexes
To identify sex-specific gene expression and genomic loci, we produced 10-20 M reads of RNA-Seq data for male and female specimens (Table S3) that mapped approximately 80-90% of the genome.Based on these data, we quantified and conducted differential gene expression analysis.Principal Component Analysis (PCA) of the expression profiles indicated a clear distinction between the male and female samples (Fig. 1A).A total of 9,015 transcripts were differentially expressed, with 4,685 and 4,329 transcripts showing higher expression in females and males, respectively (Fig. 1B).Gene ontology enrichment analysis of each gene set indicated enrichment of various pathways (Figure S1).For females, we observed enrichment of RNA processing, cellular component biogenesis, and negative regulation of biological processes.In contrast, terms related to cyclic nucleotide biosynthetic processes, aminoglycan metabolic processes, and monatomic ion transport were enriched in males.
The sex determination cascade comprises multiple genes, forming a signaling cascade that causes differentiation between the sexes.We first focused on DMRT, a well-conserved gene family that regulates sex-specific tissue development and behavior.Initial BLAST analysis identified five DMRT orthologs (PARRI_0009851, PARRI_0001169, PARRI_0005877, PARRI_0003090, and PARRI_0003093) [18].We observed that three genes, PARRI_0003090, PARRI_0003090, and PARRI_0005877, were upregulated in males, whereas PARRI_0003090 was moderately expressed (TPM > 30) in males.
A diverse array of lineage-specific upstream signaling factors (e.g., tra2, nix, fem) induce sex-specific splicing variants of the doublesex gene, transmitting signals to the downstream sex development cascade [6].Although the master regulator of sex determination is highly variable, several components of the cascade are highly conserved, such as the DMRT orthologs.Likewise, the transformer-2 (tra2) gene is a DNA-binding protein coupled with the Tra protein, causing sex-specific splicing of dsx in insects [6].The P. metropolitanus Tra2 (PmTra2; PARRI_0000692) exhibited sex-biased variants (Figure S2).We also sought to identify genes participating in the sex cascade in Drosophila, e.g., fruitless (fru) and sexlethal (sxl) genes.BLAST searches identified two candidates for sxl orthologs (PARRI_0002227 and PARRI_0007430) [18], which showed contrasting expression profiles.However, a phylogenetic analysis of these genes could not determine whether these orthologs were sxl or a gene family with relatively high similarity; therefore, we cannot conclude whether these genes are sxl orthologs (data not shown).No hits were found for fru in P. metropolitanus nor any tardigrade genomes.Thus, we concluded that fru is missing and sxl remains questionable in P. metropolitanus.

Genomic loci of the male-biased genes
We detected a peculiar population of genes that were expressed approximately > 25 higher in males (Fig. 1B).Hypothesizing that these male-biased expressed transcripts may be sex-specific genes located on the sex chromosome, we conducted a genomic enrichment analysis to determine genomic loci enriched in these highly biased genes.Using a genomic bin of 200 kbp (corresponding to roughly 30 genes per bin) against differentially expressed transcripts that had over 10 × fold change than the other sex (Female: 674, Male: 1724), we detected 325 (29 scaffolds) and 12 (three scaffolds) bins for males and females, respectively (Fig. 2A).We noted that approximately 2% of the male-biased genes had more than TPM 10 in females (11% for male-expression of female-biased genes), thus implying the specificity of male-biased genes.Gene ontology enrichment analysis of genes located in these bins indicated a high enrichment of transcripts related to sperm function (Table S4, S5).Interestingly, two of the three male-induced DMRT paralogs (PARRI_0003090 and PARRI_0003093) were located within a bin enriched for male-biased genes on the scaffold Parri_scaffold0000005 (Fig. 2B) [18].We also observed that the genes within and in the surrounding regions of these bins were also expressed in females, suggesting that these genomic loci may not be male-specific (Fig. 2C).To evaluate whether these genomic loci enriched for male-biased genes were on the same chromosome, we preformed synteny analysis with the recently reported chromosome-level genome assembly of H. exemplaris.Paramacrobiotus metropolitanus has been observed to have a 2n = 10 karyotype, similar to that of H. exemplaris [23,75].While the queried 29 male-biased scaffolds did not focus on a particular chromosome, we observed a slight bias toward chromosomes 1, 2, and 5 (Fig. 2D).Furthermore, we observed that the region harboring the paralogous DMRT loci and the surrounding region on Parri_scaffold0000005 were missing in H. exemplaris, where a genomic region on a different scaffold (Parri_ scaffold0000002) was inserted into in H. exemplaris (Fig. 2E).These data suggest that this genomic region may have emerged in the P. metropolitanus lineage.

Emergence of a novel dmrt93B-like subfamily specific to Macrobiotidae
Given the importance of paralogous DMRT genes located on Parri_scaffold0000005, we focused on the characterization of the orthologs to determine the characteristics of these paralogs.
First, we submitted the amino acid sequences of the P. metropolitanus DMRT family for phylogenetic analysis, incorporating various tardigrade DMRT orthologs from genome and transcriptome assemblies.Careful examination of the PARRI_0003093 gene structure revealed a misassembly of a single nucleotide insertion, This caused a frameshift in the 3′ terminus, leading to a truncated coding sequence.Therefore, manual curation for this gene was performed, resulting in 463 amino acid sequences.Phylogenetic analysis identified PmDm-rt99B (PARRI_0009851), PmDmrt93B (PARRI_0005877), and PmDmrt11E (PARRI_0001169) orthologs, as well as two Dmrt93B paralogs (PmDmrt3090 PARRI_0003090; PmDmrt3093 PARRI_0003093; the 3090/3093 complex).The 3090/3093 complex contained DMRT genes only from Macrobiotidae species, suggesting the acquisition of this subfamily in this lineage (Fig. 3A, Table S1).We also observed a phylum-wide loss of the doublesex subfamily.Furthermore, we observed an Echiniscidaespecific Dmrt93B subfamily that was not included in the 3090/3093 complex, while the relatively lower bootstrap  S1).We only found the DM domain, but not the CUE-DMA domain in these subfamily members through InterProScan analysis.Interestingly, phylogenetic analysis indicated that the dsxlike gene of the velvet worm branched into a Doublesex clade with arthropods, suggesting that dsx emerged after the divergence of Tardigrada (Fig. 3A, Table S1).We did not detect two copies of the 3090/3093 complex in several other Macrobiotidae species.A direct comparison between PmDmrt3090 and PmDmrt3093 amino acid sequences indicated that the first 30-180 aa sequences were extremely similar, but the intron nucleotide sequences were completely different (Figure S4AB).Furthermore, multiple nanopore reads spanned the entire length of each gene.Together, we suggest that the two copies were not the result of misassembly of these loci.We also noted that no ONT reads spanned both PmD-mrt3090 and PmDmrt3093.However, the 3090/3093 complex region spanned more than 30 kbp and the N50 length of the ONT data was approximately 17 kbp.It is possible that there are no ONT reads spanning the entire region.While reassembly of the ONT reads using more recent assembly methods produced a more contiguous assembly (NextDenovo + NextPolish; Table S6), these two genes were predicted to be two separate genes.
Based on these annotations, we identified PmD-mrt3090, PmDmrt3093, and PmDmrt93B to be significantly expressed in males; thus, all three induced copies belong to the Dmrt93B clade (Fig. 3B).We then utilized our previously reported single specimen RNA-Seq data of the embryonic and juvenile life stages of the parthenogenetic tardigrades H. exemplaris and R. varieornatus to evaluate whether DMRT orthologs in others may be functional.Only females have been observed in both species, suggesting the lack of masculinization in these species.All three Dmrt11E, Dmrt93B, and Dmrt99B orthologs in H. exemplaris and R. varieornatus (RvDmrt11E: g5527, RvDmrt93B: g9000, RvDmrt99B: g7078; HeDmrt11E: BV898_08851, HeDmrt93B: BV898_13063, HeDmrt99B: BV898_01934.)were expressed during embryonic stages (Figure S3), where Dmrt11E preceded Dmrt99B in both species, and the three DMRT genes were expressed at lower levels in juvenile and adult stages.
We further investigated the functionality of the DMRT orthologs by functional domain detection (Fig. 3C).While all five DMRT copies harbored the DM domain at the N-terminus, they did not contain the dimerization domain known to exist in dsx proteins required for DNA binding and sex-specific splicing variants.Furthermore, we did not find the ubiquitin bindingrelated CUE-DMA domain in PmDmrt3090 by domain search analysis.However, multiple alignment of the five orthologs and the D. melanogaster DMRT sequence suggested the conservation of several residues within the region corresponding to the CUE-DMA domain, implying the conservation of this domain.By modeling the protein structure with AlphaFold2 and aligning the D. melanogaster Dmrt93B structure, we observed that the C-terminal region showed structural homology with CUE-DMA domain-like helices (RMSD: 0.276-0.574,Figure S4CDE), suggesting that PmDmrt3090 may also harbor the CUE-DMA domain.

Contradictory data from whole genome sequencing and PCR-based genotyping
Based on our observations of several male-biased but not male-specific genomic regions, we hypothesized that these regions were not sex-specific chromosome structures.To evaluate this, we sequenced the genomes of both sexes at low coverage.We produced approximately 50-60 M reads, corresponding to roughly 20-25 × coverage (Table S2).Approximately 80-90% of the reads were mapped to the genome, resulting in roughly 15-20 × coverage.
We first calculated the coverage of the 10 kbp bins genome wide.Initial PCA of the coverage profiles did not show a clear difference between males and females (Fig. 4A).We identified several bins with half of the average genome-wide coverage that were not found in females (Fig. 4B).These characteristics are similar to those of heterozygotic chromosomes, particularly the X chromosome of males in the XY sex determination system.All of the bins that were identified as male-biased by the transcriptome analysis had genome-wide average coverage, suggesting that all regions exist in females (Fig. 4C).We also evaluated whether we could detect male or female specific regions through k-mer based analysis using the YGS method.Scaffolds that have a high number of "percent validated single-copy unmatched k-mers" (P_VSC_UK) indicate sex-specificity.Although no scaffolds had a P_VSC_UK ratio of 100, we detected five scaffolds fitting the XY sex chromosome structure rather than the ZW scheme with an arbitrary threshold of P_VSC_UK > 80 (Fig. 4D, E, F).Similar profiles were observed by SPADES reassembly using all gDNA-Seq data from our and previous studies (Fig. 4G, H, I).We also noticed that many scaffolds from both assemblies had P_VSC_UK values of approximately 50% in both female-to-male (XY) and male-to-female (ZW) analyses (Fig. 4F, I), which indicates that the corresponding region is both male-and female-specific.
To evaluate the male specificity observed in the in-silico analysis, we designed several primers to amplify regions in the scaffold Parri_scaffold0000295 that were identified as male-specific (Fig. 4J, Table S2).Evaluating individual genomic coverage indicated that in this scaffold, a single male sample had near-zero coverage, in contradiction to the other two male samples (Fig. 4J).However, PCR genotyping indicated the existence of this region in females as well, which contradicted the results obtained from the in-silico analysis (Fig. 4K).Thus, we concluded that we could not derive sex chromosomes or malespecific regions, and the male-specific regions detected above may have been an artifact of differences between individuals.We also evaluated the male-specificity of the paralogous DMRT loci on Parri_scaffold0000005, where coverage analysis suggested that this region was not male-specific (Fig. 4L).
We noticed a relatively low level of RNA-Seq mappability to the reported genome (~ 90%), which led us to re-evaluate the current genome.Completeness analysis indicated 72.9% completeness with the most recent Fig. 4 Paramacrobiotus metropolitanus lacks sex specific regions for both sexes.A PCA of genomic coverage profiles for male and female gDNA-Seq data.B Average coverage for 10 kbp bins, normalized by the median of all bins.The brown and red lines indicate the whole-genome average and half-genome average for males and females, respectively.The range between the 0 × and 10 × coverage ratio to the median is shown.C Genome coverage of male-biased bins.D, E, G, H YGS analysis for (D, G) female-to-male (XY system) and (E, H) male-to-female (ZW system) based on the (D, E) published P. metropolitanus assembly and (G, H) SPADES reassembly.The red line indicates a P-VSC-UK threshold of 80. F, I Scatter plot for P-VSC-UK values for male-to-female and female-to-male analysis for the (F) published genome and the (I) SPADES reassembly.Contigs shorter than 1,000 bp were removed from the SPADES plot.J gDNA-Seq coverage for all samples and the location of genotyping primers within contig Parri_scaffold0000295.Blue and red correspond to male and female samples, respectively.K Electrophoresis of genotyping primers designed for male specificity.Male specificity was not observed for any primer set.L gDNA-Seq coverage of PmDmrt3090/PmDmrt3093 harboring scaffold Parri_scaffold0000005 for all samples.Colors indicate samples for each sex version of BUSCO.Furthermore, we observed several scaffolds with inconsistent coverage distribution in our sex-separated data, but not in Hara et al.Illumina data [18].Therefore, we tested if recent assemblers would result in a more contiguous and complete assembly, compared to the Canu assembled current genome.However, we were not able to obtain a more complete genome, with the maximum being a 0.4% increase for the assembly derived with NextDenovo + NextPolish (Table S6).Other statistics had a large increase; N50 from 1.0 M to 1.3 M, longest scaffold length 4.48 M to 9.23 M. For comparison, we evaluated the completeness of other highquality tardigrades genomes, namely R. varieornatus and H. exemplaris.Both BUSCO and compleasm resulted in completeness values similar to P. metropolitanus, R. varieornatus (C:74.6%)and H. exemplaris (C:73.3%).These data suggests either tardigrade genomes may lack some BUSCO genes, or the gene detection algorithm of the current BUSCO software may not fit the genome of tardigrades, resulting in lower BUSCO scores.Therefore, we used the current genome for P. metropolitanus for later analysis.

Sex-bias in anhydrobiosis-related genes
A major feature of tardigrades is their ability to survive environmental extremes, a phenomenon known as cryptobiosis [76].Tolerance to near-complete desiccation is known as anhydrobiosis [77].Several tardigrade-specific gene families, i.e. cytosolic-abundant heat soluble (CAHS) and Secretory Abundant Heat Soluble (SAHS), have been implicated in anhydrobiosis protection [22].A recent study observed tissue-specific expression of anhydrobiosis genes [78].Both males and females are capable of anhydrobiosis, in which protective genes are expressed in sex-specific organs, such as the testes or ovaries.Therefore, we hypothesized the presence of sex-biased anhydrobiosis genes.
We used our previously reported RNA-Seq data for the hydrated active state and the tun state, desiccated for two days, to identify genes induced during anhydrobiosis.We detected approximately 4,500 differentially expressed transcripts, slightly fewer than in our previous report, possibly due to the different methods used for differential expression analysis.We then compared the expression profiles of anhydrobiosis and between sexes and observed approximately 1,800 transcripts that were differentially expressed under both conditions (Fig. 5A).As hypothesized, we observed that three CAHS and one SAHS ortholog were sex-biased, possibly indicating tissue specificity (Fig. 5A).Interestingly, all three CAHS orthologs induced in males were the only three among the 13 CAHS orthologs that were not differentially expressed during anhydrobiosis (Table S7).Phylogenetic analysis indicated that these CAHS orthologs were CAHS1 (PARRI_0016931), putative CAHS5 (PARRI_0006576), and CAHS5 (PARRI_0002229) orthologs, following the proposed naming scheme of Fleming et al. [54].In contrast, the SAHS ortholog, detected as differentially expressed, was induced in the females.We also found six orthologs of tardigrade-specific manganese-dependent peroxidase [33] to be highly expressed in males but not in females.Only four genes were found to be induced during anhydrobiosis.

Discussion
In this study, we focused on gonochoristic tardigrade P. metropolitanus to identify possible factors that affect sexual dimorphism.Cytological studies have not identified definitive sex-linked chromosomes in tardigrades [80,81], and multiple reports have observed biased sex ratios in tardigrades [23,[82][83][84][85].These observations suggests that sex determination in tardigrades may not depend on the random distribution of sex chromosomes (or the existence of a sex chromosome).Even in the absence of sex chromosomes, as hypothesized in tardigrades, genomic loci affecting sexual dimorphism could exist, which may be detected by comprehensive omics methods.
Therefore, we aimed to characterize the molecular basis of sexual dimorphism in tardigrades by comparing the transcriptome and the genome between the sexes of P. metropolitanus.We hypothesized that sexlinked genes may be related to sex determination or dimorphism, and if focused on a small genomic region, may imply a sex-determining region, such as the M factor found in many eukaryotes [86].Transcriptome analysis of both sexes indicated a large number of sex-biased genes, despite the small morphological sex-linked differences in Macrobiotidae, with the exception of their germline [20].We observed upregulation of genes related to spermatogenesis in males, which reflects the activation of spermatogenesis, and large amounts of sperm are continuously produced in adult males [23,83].In contrast to that in males, DNA replication-and meiosisrelated genes were highly expressed in females.Females undergo DNA replication not only to produce oocytes through meiosis [23] but also to shed the cuticular exoskeleton during the last stage of the reproductive process (simplex stage) [23,87].Mitotic cells are generally observed in the post-simplex stage [88].Together, the regulation of DNA replication and meiosis is consistent with the production of mitotic cells and extensive replication of the epidermal layer [88,89].
We identified a small gene set highly biased toward males, but missing in females, which we hypothesized may be related to sexual dimorphism.Genome loci enrichment analysis of this gene set identified approximately 325 bins spanning 29 scaffolds as male-biased.This region was enriched in sperm and ion transportrelated genes, which is consistent with the production of sperm at the adult male life stage.To evaluate sex specificity, we produced low-coverage genome sequencing data to evaluate sex-specific regions and observed that most regions were present in the genomes of both sexes.Genome-wide analysis revealed several male-specific regions; however, PCR evaluation produced contradictory results.We used a laboratory-cultured TYO strain of P. metropolitanus for genome and transcriptome sequencing, therefore, we anticipated low levels of heterozygosity within the culture population.However, the results obtained at this stage suggest that the genomic differences we detected as sex-linked can be explained by individual variability.Additionally, during the YGS analysis, we observed a high number of contigs that showed approximately 50% P_VSC_UK, suggesting that there are a large number of contigs that contain sequences specific for both sexes, which we hypothesize that individual variability may have caused this abnormal distribution.Together, the lack of sex-specific regions may indicate that the difference between sexes is due to epigenetic modifications.
One of the key findings of this study is the accumulation of knowledge for sex determination cascade-related genes, particularly the DMRT gene family.The DMRT family is a highly conserved transcription factor that plays an important role in sex differentiation in many animals and has been studied extensively in insects [7].Several studies have identified DMRT orthologs to be located on the sex chromosomes and regulate the growth of sex-specific tissues [8,9].The evolutionary background of this gene family has been extensively analyzed in other lineages [7]; however, such analysis has been overlooked.In our analysis, we identified a Macrobiotidae-specific Dmrt93B subfamily located in a male-biased region, which we termed the 3090/3093 complex in addition to the Dmrt99E, Dmrt93B, and Drmt11E subfamilies.Orthologs of DMRT genes in both gonochoristic and parthenogenetic tardigrades are expressed during several stage of development, thus the genes may be functional.Family-specific orthologs of Dmrt93B subfamily have been found in Macrobiotidae and several Echiniscidae.The lack of two copies in other Macrobiotidae species may be the result of misassembly in their genomes; the analyzed genomes are based on Illumina short reads, and the extremely similar 30-180 aa (corresponding to approximately 450 bp) may have resulted in a misassembly.While conservation in Echiniscidae complicates the evolution of this subfamily, the identification of orthologs in various Macrobiotidae species suggests that this is an important DMRT subfamily.In fact, the two 3090/3093 complex paralogs were expressed higher in males, similar to Daphnia Dsx1 [11,12], suggesting that these subfamily orthologs may inhibit feminization or progress musculation.Furthermore, we did not find any orthologs of the dsx subfamily in any of the tardigrade genomes analyzed, and we did not identify splicing variants in any of the P. metropolitanus DMRT orthologs, suggesting a sex differentiation cascade different from those that rely on sex-specific dsx splicing variants like those observed in insects.In addition, we observed sex-biased expression in tra2 splicing variants which may be functional in the P. metropolitanus sex determination cascade; however, several factors in this cascade may be lost in this lineage.The Bombyx mori sxl gene induces dimorphism of the sperm, not sex determination [90]; therefore, it is possible that the lack of sxl may imply a different regulatory pathway than is known.
Tardigrades are renowned for their ability to tolerate extreme stress [22], and P. metropolitanus also shows a high tolerance to desiccation stress [18].Interestingly, we observed sex-biased expression of several anhydrobiosis genes, hypothesized to play protective roles during anhydrobiosis [30][31][32][33]91].For instance, CAHS genes are tardigrade-specific proteins that form gel filaments that possibly protect cells [92][93][94].Recent studies have observed tissue/organelle specificity for these proteins, which further implies the existence of orthologs with sex-specific expression [78].Therefore, we hypothesized that orthologs of such genes may exhibit sex-specific expression to protect sex-specific organs.Indeed, we identified CAHS, SAHS, and AMNP orthologs with sex-specific expression.Two of the three male-induced CAHS orthologs were highly expressed, but were not induced differentially between active and anhydrobiotic conditions.This may imply that these CAHS orthologs participate in the protection of male-specific tissues or sperm.Furthermore, we identified a P. metropolitanus Dsup ortholog that is highly expressed in females.Coupled with the observation of the enrichment of meiosisrelated genes from transcriptome analysis, we suggest that Dsup may actively function to accommodate the production of oocytes/oogenesis rather than spermatozoa/spermatogenesis.In contrast, AMNP, a tardigradespecific peroxidase, was highly expressed in males, suggesting enhanced protection against oxidative stress.Similar observations have been made in the sperm of many animals [95,96].Together, the sex-biased expression of anhydrobiosis genes may provide protection for sex-specific tissues.

Conclusions
In this study, we identified male-biased regions that may harbor potential candidates that regulate sexual dimorphism in the gonochoristic tardigrade P. metropolitanus.Simultaneously, these data denied the sex-chromosomebased sex determination scheme.We also provide evidence for a new DMRT subfamily that may contribute to sex differentiation in this family.The 3090/3093 complex DMRT paralogs may be initial candidates for disruption or gene editing for evaluation their relationships with sex determination [78,[97][98][99].Future studies utilizing highquality genomes and careful physiological experiments are required to reveal sex determination cues not only in this species but also in other tardigrades.

Fig. 1
Fig. 1 Transcriptomic analysis of both sex.A PCA of expression profiles.B Scatterplot of the expression profiles.Red dots indicate differentially expressed transcripts (FDR < 0.05)

Fig. 2
Fig.2Multiple male-biased regions within the P. metropolitanus genome and their synteny.A Genome-wide enrichment analysis of male-or female-biased transcripts.Scaffolds were ordered by size and colored green and purple to visualize the scaffolds.The threshold of FDR < 0.01 was used.B, C Characteristics of scaffold Parri_scaffold0000005 harboring the DMRT paralogs.B FDR values from the genomic loci enrichment analysis plotted against the bin's position.Blue and red indicate the values for males and females, respectively.C Expression fold-change log2(male + 0.1 / female + 0.1) of the genes plotted against their locations along the scaffold.Colors indicate whether the gene was differentially expressed.D Macro-scale synteny analysis to identify orthologous genomic loci in male-biased scaffolds.Gray lines indicate syntenic blocks between H. exemplaris and P. metropolitanus scaffolds.The synteny block highlighted in green indicates the location of DMRT paralog loci.The numbers on the bar indicate the chromosome number or scaffold ID for each genome assembly.E Synteny region of the DMRT paralog loci on P. metropolitanus scaffold Parri_scaffold0000005 and H. exemplaris Chromosome 1. Orthologs between H. exemplaris and P. metropolitanus determined by synteny analysis are connected by a gray ribbon.The corresponding region in H. exemplaris Chromosome 1 also aligned to the P. metropolitanus scaffold scaffold0000002 3.79-3.84Mb, as indicated in between He Chr1 10.02-10.23 Mb and Pm scaffold000005 1.96-2.37Mb.The Dmrt3090/3093 complex orthologs (PmDmrt3090 and PmDmrt3093) are also indicated

Fig. 3
Fig. 3 Phylogenetic analysis and the expression of DMRT orthologs.A Phylogenetic analysis of DMRT orthologs detected in the tardigrade genomes.The DMRT families were classified based on the orthologs of the model species.Bootstrap values of less than 90% are shown on the branch.The blue colored region indicates the Macrobiotidae-specific orthologs.The P. metropolitanus and Echiniscidae-specific Dmrt93B orthologs (Cornechiniscus lobatus, Echiniscus spp.Pseudoechiniscus sp.), and the partial dsx of velvet worm (Euperipatoides kanangrensis) are indicated in bold.B Expression of PmDmrt orthologs.Triangle points indicate differentially expressed genes and circles indicate non-significant changes.The gray line indicates x = y.C Multiple alignments of DM and CUE-DMA domains.Dm indicates D. melanogaster

Fig. 5
Fig. 5 Sexual bias in anhydrobiosis genes and identification of PmDsup ortholog.A Comparison of gene expression profiles between the sexes during anhydrobiosis.Log2 (Tun + 0.1) / (Active + 0.1) were plotted for the x-axis and for the y-axis log2 (Male + 0.1) / (Female + 0.1).Red dots indicate transcripts detected as differentially expressed in both comparisons.Sex-biased, but not anhydrobiosis induced, anhydrobiosis related genes have been noted by colored annotations: Blue: CAHS, Green: SAHS, Red: AMNP, and Black: Dsup.B Synteny analysis to identify orthologous genomic loci in H. exemplaris and P. metropolitanus.HeDsup and PmDsup have been annotated in the plot.C DISOPRED and IUPRED3 scores (D) Protein structure predicted by ColabFold.The N-terminus to the C-terminus shows gradient colors from blue to red.E Expression of PmDsup in both sexes